Diffusion on random site percolation clusters. Theory and NMR microscopy 

experiments with model objects 
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Quasi two-dimensional random site percolation model objects were fabricated based on computer 
generated templates. Samples consisting of two compartments, a reservoir of H2O gel attached 
to a percolation model object which was initially filled with D2O, were examined with NMR (nu- 
clear magnetic resonance) microscopy for rendering proton spin density maps. The propagating 
proton/deuteron inter-diffusion profiles were recorded and evaluated with respect to anomalous dif- 
fusion parameters. The deviation of the concentration profiles from those expected for unobstructed 
diffusion directly reflects the anomaly of the propagator for diffusion on a percolation cluster. The 
fractal dimension of the random walk, d w , evaluated from the diffusion measurements on the one 
hand and the fractal dimension, df, deduced from the spin density map of the percolation object on 
the other permits one to experimentally compare dynamical and static exponents. Approximate cal- 
culations of the propagator are given on the basis of the fractional diffusion equation. Furthermore, 
the ordinary diffusion equation was solved numerically for the corresponding initial and bound- 
ary conditions for comparison. The anomalous diffusion constant was evaluated and is compared 
to the Brownian case. Some ad hoc correction of the propagator is shown to pay tribute to the 
finiteness of the system. In this way, anomalous solutions of the fractional diffusion equation could 
experimentally be verified for the first time. 

PACS numbers: 05.40.-a, 82.56.Lz, 47.53.+n, 64.60. Ht 



I. INTRODUCTION 

Randomly disordered media are present in many fields 
of nature and science. The dynamical properties ruled by 
the geometrical structure are of special interest in fields 
of physical and engineering processes like filtering and 
exploration jj], ||, [| . Percolation theory has proven to be 
a powerful tool to model porous systems || [l6|, [DJ . 

The objective of this study is to examine diffusion on 
random site percolation clusters experimentally and an- 
alytically. There are several numerical simulation stud- 
ies in the literature suggesting an anomalous displace- 
ment behavior related to the fractal nature of the clusters 
H ||, 0]. However, there is little experimental evidence 
for the reality and practical detectability of anomalous 
diffusion so far |,§ @. 

The objective of the present work is to exploit a new 
experimental strategy. This is: (a) to generate numeri- 
cally a percolation model cluster, (b) to determine the 
characteristic parameters numerically, (c) to fabricate 
model objects using the percolation clusters as templates, 
(d) to record nuclear magnetic resonance (NMR) spin 
density maps from the (water-filled) pore space, (e) to 
evaluate the characteristic cluster parameters on this ba- 
sis again, (f) to study intcrdiffusion of heavy and light 
water in the pore space, and (g) to compare the exper- 
imental interdiffusion profiles with solutions of the frac- 
tional diffusion equation |rT[| for the first time. In a sense, 
we are thus continuing our previous work in which we had 
already explored static and dynamic properties in various 
three-dimensional and quasi two-dimensional percolation 
model objects @ [l| [b§ fjj). 



Random site percolation structures are defined in the 
two-dimensional case by sites on a square lattice. They 
are occupied with a probability p that is usually cho- 
sen in the vicinity of the (two-dimensional) percolation 
threshold, p c = 0.592746. Neighboring occupied sites are 
connected by pores with a cross-section corresponding to 
the lattice constant or integer multiples of it. The total 
subset of connected lattice sites form a so-called clus- 
ter. For p>p c , sample-spanning clusters occur that can 
be examined with respect to transport properties. The 
pore space structure generated by the random site per- 
colation model can be characterized by four parameters, 
that is, the lattice constant, a, the fractal dimension, 
df, the correlation l eng th, £, and the percolation prob- 
ability, Poo Jl6|, [I?], The latter quantity is defined 
as the probability that a site belongs to the "infinite" 
cluster traversing the whole sample |19| . The correlation 
length, which is of particular interest here, is defined as 
the mean distance between two sites of a finite cluster 
(or the mean hole diameter in an infinite cluster). In the 
real percolation model objects we are considering here, 
the minimum lattice constant (or pore diameter) is given 
by the mechanical resolution of the fabrication process 
(see below). 

Random site percolation clusters are known to display 
fractal properties on a length scale below the correlation 
length. That is, the volume-averaged porosity scales with 
the probe volume radius, r„, as 



/ E for a < r p < f (l) 
P 1 "ooocro for r p > f. W 

The Euclidean dimension is denoted by d,E (= 2 in the 
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present case). The fractal dimension for 4 = 2 was 
theoretically derived as d f = 91/48 ps 1.896 |§. 

The purely structural relationship Eq. ([!]) is opposed 
by the dynamic property for the mean squared displace- 
ment of a random walker on the cluster, 



t 2/d w 
Defft 



for 
for 



t a -C t < t£ 
t > t(, 



(2) 



where oc £ d ™ is the time the random walker needs to 
explore the correlation length £, and d w is the fractal di- 
mension of the random walk. The lower time limit of the 
anomalous diffusion regime is given by the time needed 
for the displacement of length a, t a . The diffusion coef- 
ficient becoming effective in the long-time limit, t 3> t^, 
is denoted by D e ff. According to the Alexander/ Orbach 
conjecture (2C|, the quantity d w is assumed to be related 
to the fractal dimension as 



d w = -dt 
2 



for 



d E > 2. 



(3) 



That is, the structural parameter df characteristic for 
the volume-averaged porosity is linked to the dynamic 
parameter d w specifying anomalous diffusion. For c2e = 
2, the diffusion exponent becomes d w ~ 2.87. In this 
study experimental evaluations for both quantities have 
been carried out, so that a comparison becomes possible. 
Note, however, that Eq. (ft) is not considered to be an 
exact relation Q (L7|, JyJ |§ . 

A theoretical problem of intriguing impact is the 
complete propagator description of anomalous diffusion 
rather than restricting oneself to the second moment of 
the propagator according to Eq. (^) . In the second part of 
this article, the analytical treatment based on the frac- 
tional diffusion equation fll] , |23| is outlined and com- 
pared with the inter-diffusion profile data acquired in our 
experiments. 



II. TECHNIQUES AND INSTRUMENTS 
A. Methods for measuring diffusion 

In the model objects to be studied here, the minimum 
pore diameter is Ar = 400 //m. The displacement length 
scale needed to probe anomalous diffusion, i.e. displace- 
ments obstructed by the matrix, is r<j 3> Ar. The or- 
dinary pulsed gradient spin echo technique (see Refs. 

p5| , for instance) is therefore not suitable for the 
detection of anomalies in liquids in the present situation. 

The much larger displacement rate in gaseous phases 
would permit such studies in principle. In Ref. |ll| we 
have studied diffusion of methane gas in a percolation 
model object. Although there was some indication of 
an anomalous behavior, the experiment turned out to be 
difficult because of the poor detection sensitivity. In this 
respect, diffusion studies using laser-polarized or ther- 
mally polarized 129 Xe are more promising pi, |2^. Also, 



the use of inert fluorinated gases possibly at somewhat 
elevated pressures may be more favorable [^8[ |2!|. In 
any case, there is a diffusion mechanism ("Knudsen dif- 
fusion") relevant in gases which is different by nature 
from diffusion in liquids fl. 

We therefore preferred to employ an isotope interdif- 
fusion method. The samples consisted of two compart- 
ments initially filled with H2O (in gel form) and D2O. At 
the beginning of the experiment the compartments were 
pressed on each other in close contact, so that interdiffu- 
sion was initiated. The time evolution of the proton spin 
density maps in the D2O compartment was then stud- 
ied as a function of time. The technique was ordinary 
one- or two-dimensional NMR imaging of concentration 
profiles, an application for diffusometry purposes already 
described in Refs. 69, |3ll 51. 



The proton density profiles were recorded either in the 
form of one-dimensional spin density maps, or were eval- 
uated from the two-dimensional spin density maps by 
projection on the main diffusion direction. The latter 
variant has the advantage that signal noise from matrix 
areas can be screened off before evaluating the profile 
data. 

The propagation of the proton density profile at half 
height as a function of time permits one to determine the 
time dependence of the mean squared displacements. Al- 
ternatively, the profiles themselves at a given time can be 
examined with respect to the character of the diffusion 
process. In the latter case, the full propagator character- 
istics and not just its second moment are mattering. 



B. NMR tomograph and acquisition parameters 

The one- or two-dimensional proton density maps of 
the water filled pore space of percolation model objects 
were recorded with the aid of a NMR tomograph con- 
sisting of a 4.7 T Bruker magnet with 40 cm horizon- 
tal room temperature bore and a home made radio fre- 
quency console. Typical radio frequency and field gra- 
dient pulse schemes for spin-echo NMR imaging can be 
found in Ref. |24j, for instance. The spatial resolution of 
the images was better than 300 /im. The acquisition of 
a two-dimensional spin density map typically took 20' to 
60', so that a reasonable time resolution was given. 

Isotopic dilution by deuterons prolongs the local trans- 
verse and longitudinal relaxation times due to the re- 
duced number of dipolar interaction partners (see Ref. 
p3| , for instance). For the evaluation of spin density 
maps, the spin echo signals therefore have to be corrected 
if the repetition time is not much longer than the longest 
proton spin-lattice relaxation time, Ti, or if the echo time 
is not much shorter than the transverse relaxation time, 
T 2 . 

Typical echo times, Tg, were between 20 and 30 ms. 
This is to be compared with transverse relaxation times 
of several seconds in water at room temperature. Signal 
attenuation on this basis is therefore totally negligible. 
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The situation is less clear with the effect of spin-lattice 
relaxation. The repetition time, Tr, typically was 2 s, so 
that the spin density profiles could be distorted at the 
low-concentration side by saturation effects. In some of 
the experiments, we have therefore varied the repetition 
time between 0.25 and 12.2 s in order to evaluate the local 
spin-lattice relaxation times. The local signal intensities 
were then corrected correspondingly to provide the true 
spin density profiles. No significant spin-lattice relax- 
ation effect could be diagnosed (see the data discussed 
below) . 

C. Computer generated percolation clusters 

In the insets of Fig. || and Fig. |], typical two- 
dimensional random site percolation clusters generated 
on a square lattice are shown. The occupation probabil- 
ity, p, is slightly above the percolation threshold value 
p c = 0.5927 for the Euclidean dimension cIe = 2 |m[ . 
The volume-averaged porosity was evaluated using the 
so-called sandbox method (l^, |l3): N p probe circles of 
varying radius r p are first placed randomly at positions 
rk within the cluster in such a way that the center of 
the probe volume (which actually is an area in the two- 
dimensional case) is in the pore space. Then the average 
values of the observables are formed for the Nv voxels at 
positions rj inside the probe volume. Finally, the arith- 
metic mean of the data set for the N p probe volumes with 
a given radius r p is taken. In other words, the volume- 
averaged porosity is defined as 

1 N p N v 

p k=i v j=i 

where r > |r& — rj\. This quantity can also be evalu- 
ated from black-and-white converted, experimental spin 
density maps as described in Ref. |lj| . 



matrix is unfavorable because of the tendency of forming 
voids upon gelation. It turned out that in this case the 
stabilizing effect of the solid matrix is sufficient. 

The influence of deuteration and gelation on the bulk 
water self-diffusion coefficient was checked in an ordinary 
pulsed-gradient spin echo experiment p4]. At 20 °C and 
a gel content of 1.2 %, the self-diffusion coefficients of 
H 2 and D 2 were found to be 1.8 x 10 -9 m 2 /s and 1.4 x 
10~ 9 m 2 /s, respectively. That is, the water self-diffusion 
coefficient is slightly reduced both by deuteration and 
gelation. 

An isotope effect is known to show up already without 
gelation: Mills [[35] reports self diffusion coefficients in 
pure H 2 and D 2 at 25 °C D H . l0 = 2.3 x 10~ 9 m 2 /s 
and Du 2 o = 1-87 x 10~ 9 m 2 /s, respectively. That is, 
smaller diffusion coefficients are expected with increasing 
dilution of 1 H as it occurs at the diffusion front. 

Furthermore, there may be a difference in the chemi- 
cal potential in pure H 2 and pure D 2 0, so that a slight 
deviation from the self-diffusion situation may play a 
role. On the other hand, no significant influence on the 
shape and time evolution of the inter-diffusion profiles 
was found in the bulk-to-bulk experiments described be- 
low. The conclusion therefore is all these effects are of 
minor importance or largely compensate each other so 
that (partial) deuteration and gelation does not percep- 
tibly affect the percolation cluster characteristics of the 
propagator to be probed. 

The quasi two-dimensional model objects were kept in 
a horizontal position during the whole measuring process 
in order to avoid any convective displacements. In or- 
der to improve the detection sensitivity, several identical 
slices were stacked on each other (see Fig. [I]) . All experi- 
ments were carried out at room temperature, (21± 1) °C. 

III. RESULTS 
A. Volume-averaged porosity 



D. Model objects and measuring conditions 



The percolation model objects were fabricated using a 



circuit board plotter (for details see Refs. J_2| 
based on the computer-generated templates. 
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The me- 
chanical fabrication resolution was Ar = 400 /im (see 
photograph in Fig. |l|). The adjusted milling depth ranged 
from 1 to 2 mm in the different objects produced. 

The objects were filled with heavy water and brought 
into contact with reservoirs of H 2 gel (Kelcogel, 1.5 % 
by weight) at time t = 0, when the interdiffusion process 
was to begin. The reservoirs are schematically shown in 
the insets of Fig. ^ and Fig. ||. The proton distribution 
in the objects was then measured as a function of time 
in the form of spin density maps as described above. 

The gel form of the undeuterated moiety of the sample 
was needed for stabilization and in order to prevent flow. 
On the other hand, gel stabilization inside the percolation 



The volume-averaged porosity was evaluated as a func- 
tion of the probe volume radius on the basis of Eq 
for the percolation clusters shown in the insets of Fig, 
and Fig. |sj. Corresponding evaluations from experimen- 
tal spin density maps lead to equivalent decays as demon- 
strated in Ref. Jig ]. The fractal dimension according to 
Eq. (|J) was found to be df = 1.87 in both cases. This 
value will be compared with the experimental value for 
d w according to Eq. (0). 



B. Isotope concentration profiles for interdiffusion 
between bulk water compartments 

In order to test our measuring and evaluation tech- 
nique, we have recorded isotopic inter-diffusion profiles 
in a sample consisting of two equal compartments ini- 
tially filled with bulk H 2 and D 2 gels, respectively. In 
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this case, the initial distribution of the proton density is 
of the type 



C(x,t = 0) 



C 




for 
for 



x < 
x > 



(5) 



Provided that the diffusion process is normal, that is, in 
the absence of obstructions by a matrix, thejaroton spin 
density profiles at later times are given by 



f 



C(x,i) = -Cq erfc 



1-JDt 



(6) 



where D is the diffusion coefficient and erfc(z) is the com- 
plementary error function. In principle, this solution of 
the ordinary diffusion equation applies to "infinite" sys- 
tems, that is, to root mean squared displacements much 
less than the extension of the sample. 

Figure ^| shows experimental proton spin density maps 
measured in a two-compartment sample (for an illustra- 
tion see inset of Fig. ||) as a function of time t after bring- 
ing the compartments into contact with each other. The 
right and left compartments were initially filled with bulk 
H 2 and D 2 gels, respectively. The mean proton spin 
density profiles, i.e. projections of the two-dimensional 
spin density maps on the main diffusion coordinate axis, 
x, are also shown in diagram form (white lines). 

In Fig. ||, these experimental concentration profiles are 
compared with those predicted by Eq. @. The fits of 
Eq. (^) to the experimental data reproduces the room 
temperature value of the diffusion coefficient measured in 
bulk water with the pulsed gradient spin echo technique 
fjM Ijjfj], D = 2 x 10~ 9 m 2 /s, very well with the exception 
of the shortest and longest diffusion intervals. 

At the shortest diffusion interval, imperfections of the 
initial isotope distribution and the limited time resolu- 
tion of the (two-dimensional imaging process) are ex- 
pected to matter. The concentration profiles at the 
longest diffusion times are already affected by the finite 
extension of the sample which conflicts with the assump- 
tion in Eq. (^) of infinite compartments. 

Taking together the potential sources of systematic ex- 
perimental errors mentioned before, one can state that 
the agreement between the theoretical and experimental 
concentration profiles in bulk samples is very reasonable 
so that reliable evaluations of anomalous diffusion fea- 
tures in the percolation model objects can be expected. 
This is corroborated by the determination of the mean 
square displacement of the diffusion front in a second ex- 
periment. The set-up is schematically shown in the inset 
of Fig. H In this case, the half-height positions of the con- 
centration profiles (projections of the two-dimensional 
spin density maps on the main diffusion direction) were 
evaluated, squared and plotted versus time. 

The middle section of the experimental curve shown in 
Fig. [|can nicely be described by a power law x\, 2 cx t 105 
which is very close to the linear mean squared displace- 
ment law for normal diffusion. As in the experiment dis- 
cussed before, the deviations at short times reflect the 



initial situation that can only imperfectly be described 
by a step function. The plateau reached at long times is 
due to the finite extension of the sample. 



C. Isotope concentration profiles for interdiffusion 
between bulk water and water filled percolation 
cluster compartments 

Equation (^|) is based on a Gaussian propagator as ex- 
pected for ordinary diffusion. That is, it does not account 
for diffusion in a percolation cluster where the second mo- 
ment of the propagator obeys an anomalous-diffusion law 
as given in Eq. (|2|). 

Figure || shows two-dimensional proton spin density 
maps acquired in an experimental set-up schematically 
shown in the inset of Fig. |[ Inter-diffusion between 
a H2O gel filled reservoir and a stack of quasi two- 
dimensional percolation model objects was examined. 
Superimposed to the spin density maps the mean spin 
density profiles (i. e. projections on the x axis) are shown 
in Fig. ^ as white lines. 

The half-height positions of these mean concentration 
profiles, xi/ 2l were evaluated, squared and plotted ver- 
sus time as shown in Fig. [| In the limit of long diffusion 
times when the imperfections of the initial isotope distri- 
bution and the finite time resolution of 26' do not mat- 
ter anymore, the data can be described by a power law 
again. The imperfection of the initial isotope distribution 
is matched by the revised origin of time t. Fitting the ex- 
ponent of Eq. (||) for t <C to the data leads to d w = 2.89 
which favorably compares to the Alexander/ Orbach con- 
jecture, Eq. (@), with angleichen the fractal dimension 
df determined from the very same percolation cluster 
according to Eq. (|l|) for a <C r p <C £. 

The isotopic interdiffusion profiles can also be obtained 
directly with the aid of one-dimensional imaging along 
the main diffusion direction. Fig. [?] renders typical pro- 
files recorded in this way from an experimental set-up 
schematically shown in the inset of Fig. ||. The data are 
corrected for saturation effects due to incomplete spin- 
lattice relaxation after the repetition interval. From these 
profiles, the mean squared proton displacement (in the 
percolation cluster moiety) can be evaluated as a func- 
tion of the diffusion time. Figure || shows corresponding 
data. Details are described in the legend. With a fitted 
exponent parameter of d w = 2.86 the anomalous diffusion 
character in the percolation cluster is again corroborated 
in good agreement with the Alexander/ Orbach conjec- 
ture Eq. (||). 



IV. ANALYTICAL PROPAGATOR 
TREATMENT BASED ON THE FRACTIONAL 
DIFFUSION EQUATION 

The propagator (or Green's function) for a Brownian 
random walk process in one dimension is, of necessity, 
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given in terms of the Gaussian 



C(x,t) = 



1 



: exp 



x 

"17 



(7) 



due to the central limit theorem pq . Expression (J7]) in- 
cludes the (5-initial condition C(x, 0) = <5(x) and fulfills 
natural boundary conditions \im^ x ^ ^ C(x,t) = 0. Wc 
use the rescaled time r = Z?i. The propagator ([?]) satis- 
fies the diffusion equation 



d 2 



(8) 



Representing the experimental set-up shown in the in- 
set of Fig. f|by the initial condition C(x, 0) = for x > 
and through the boundary condition C(0,t) = Co, i.e., 
we assume that, due to the comparatively large reservoir, 
a constant concentration is kept at the boundary x — 0. 
The solution of Eq. (||) for these initial and boundary 
conditions is given through 



C(x, r) = Coerfc 



2^ 



, x > 0. 



(9) 



The solution given in Eq. (|]) reaches the plateau 
C(x, t) = Cq for x 2 <C t. Consequently, the expression 
in Eq. (^|) is not normalized. Its mean, 

(1(t))=/ C(x,T)dx = 2C y/rfr, (10) 
Jo 

grows with the square root in time, and is proportional 
to the concentration Co- The second moment becomes 



{x 2 {r)) = a -C^' 2 l^. 



(11) 



In normalized form, that is, the concentration profiles are 
divided by the mean given in Eq. ( |To) ) , the mean squared 
displacement corresponds to 



{x 2 {r)) n = 



_ 4 
(l(r)) ~3 T - 



(12) 



It should be noted that the concentration profile, 
C(x, r), for the different realization with an infinite reser- 
voir, Eq. (^|), differs from the case with constant bound- 
ary concentration, Eq. (^J), only by a factor ^. Systems 
with constant boundary concentration and infinite reser- 
voir thus behave congruently. 

Let us now address how we can derive the analogous 
quantities for the anomalous diffusion data described in 
Section [II C. As mentioned above, the average diffusion 



in the percolation cluster close to the percolation thresh- 
old is anomalous in the sense that the mean squared dis- 
placement follows (x 2 (t)) oc T a where a = 2/d w [see 
Eq. (J|)], and d w = 2df/d s |3S[ |. Here, the rescaled time 
t = {DaY^t includes the anomalous diffusion constant 
D a of dimension [D a ] = m 2 /s a (T| g(J. Thus, the diffu- 
sion dynamics is controlled by two parameters, the frac- 
tal dimension df and the spectral dimension d s . Usually, 



the inequality d w > df is fulfilled so that the anoma- 
lous diffusion process is subdiffusive, i.e., < a < 1. In 
the experiment, the three-dimensional probability den- 
sity function is projected onto a one-dimensional proba- 
bility density function, and therefore the geometry in the 
pseudo-one-dimensional process becomes averaged. 

In the original geometry, the spreading of the random 
walker in space is slowed down in comparison to the free 
diffusion, due to the presence of bottlenecks and dead 
ends on all length scales, leading to the subdiffusive na- 
ture of the mean squared displacement. In the projec- 
tion, the trapping in a small pore of the percolation clus- 
ter corresponds to a small wiggling around some given 
coordinate. Effectively, the walker in the pseudo-one- 
dimensional measurement experiences a continued mul- 
tiple trapping process, i.e., it is immobilised for some 
"time" span r governed by the so-called waiting "time" 
distribution ip(r). According to this ip, the walker is re- 
leased and moves freely until it is trapped again. This is 
a special case of a continuous time random walk process, 
the subdiffusion being reflected in the inverse power-law 
form ip(r) ~ A a T~ x ~ a 1 42 . According to this if), the 



walker is released and moves freely until it is trapped 
again, and so forth. This stop-and-go process can be 
mapped onto the fractional diffusion equation mU E3| 



d_ 



T dx 2 ' 



C= oVl- a ^C(x,T) 



(13) 



that includes the Riemann-Liouville fractional operator 
V\- a = d/dr( V- a ) with @ 



V~ a C{x,T) 



1 



, C{x,r') 



T(a)J (t-t') 1 -"' 



(14) 



The fractional diffusion equation is equivalent to a gen- 
eralized master equation, and can be derived from a 
multiple trapping version of the fundamental Chapman- 
Kolmogorov equation |44|] . 

The Riemann-Liouville operator has the im- 
portant property that its Laplace transform is 
/ °° e~ UT T>r a f{r)dT = u- a f(u). By virtue of 
this property, it can be shown that in Laplace space, the 
solution of the fractional diffusion equation, let us call it 
C a (x, r), is connected to the solution of the Brownian 
diffusion equation (^), Cb(x,t), through the scaling 
relation C a (x,u) = u a ~ 1 CB(x,u a ) p| . This, in turn, 
can be used to find a convenient mapping between the 
Brownian solution, and any quantity obtained from it 
through linear operations^ in terms of the generalized 
Laplace transformation 



C q (x,t)= / E a (s,T)C B (x,s)ds. (15) 
Jo 



In Eq. (|15|), the function £" q (s,t) is a one-sided Levy 
distribution that can be represented in terms of Fox's H- 
function pJl]. For our purposes, we notice that E a has 
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the convenient representation 



#2/3 ( S ' T ) 



1 



T 2 /3r(i/3) 
i 



r4/3r(-l/3) 



1F1 



5 2 4s 3 \ 
6' 3'~27tV 
7 4 4s 3 



6 ' 3 ' 27t 2 



(16) 



for a = 2/3. This value for a is sufficiently close to 
the percolation cluster value 0.7, and we choose 2/3 as 
for this rational number the representation in Eq. Jl6| ) 
considerably reduces the computation time in evaluating 
the integrals of Eq. (15). In Eq. (|l6|), iF\ represents the 
confluent hypergeometric function. 

For the normalized mean squared displacement, 
Eq. (|l2|), the transformation given in Eq. ( |l5| ) yields the 
exact anomalous behaviour 



(x 2 (r)), 



T(2/3) 



2/3 



(17) 



with r = {D 2 /zf'^ 2 t [|TT|. By virtue of this expression, we 
can approximate the anomalous diffusion constant D 2 /3 
as 



D 2/3 « 1.3 x lO- 8 -^ 



(18) 



from the experimental data plotted in Fig. g where the 
observed time interval approximately fulfills the bound- 
ary condition C(0, r) — Co. To our knowledge, this is the 
first time the anomalous diffusion constant appearing in 
the above formalism has been determined experimentally. 
In the following, we use -D2/3 to construct the anomalous 
concentration profile from Eqs. (^) and (|l5|). 

In Fig. ^|, we show the anomalous concentration profile 
C a (x, t) for the longest diffusion time interval 180 h) 
in Fig. [7J, and compare it to the Brownian profile corre- 
sponding to the same interval. The latter was combined 
with the diffusion constant 2 x 10~ 9 m 2 /s (see Fig. ||). 
It is obvious that the subdiffusive profile lags far be- 
hind the wider spread of the Brownian counterpart, i.e., 
within the same time span, Brownian diffusion is more 
efficient. Quantitatively, this corresponds to the ratio 



D 2/3 / (OT(5/3)t 



l/3^| 



7 . 4t -l/3 s l/3_ 



The subdiffusive diffusion profile sequence given in 
Fig. [j] is juxtaposed to the theoretical curves in Fig. |l^ 
(a). It is obvious that the general trend follows the one 
displayed in Fig. [7|, particularly that the spacing between 
successive curves becomes less pronounced for increasing 
diffusion intervals. However, it can be realized that the 
profile falls off too fast in comparison to the experimental 
result. We believe that this is related to the fact that the 
investigated cluster is not an ideal fractal since it is finite. 
This leads to corrections in the fractional model that we 
heuristically introduce through a fudge factor as follows. 
Due to the finiteness of the systems, there exist some 
channels from the left to the right of the sample that en- 
able a much faster, essentially Brownian, exchange with 
the reservoir. 



We therefore propose an ad hoc correction, namely 
that we have two additive contributions, the anomalous 
profile plus a correction that is Brownian. That is, the 
resulting profile becomes 



C(x,t) = C a (x,T a ) + C b (x,t b ), 



(19) 



where C a corresponds to the fractional solution with 
r a = (D 2 / 3 ) 3/2 t (see Tab. | and C = 300, and C B is 
given in terms of the Brownian solution given in Eq. (j^) 
with t = Dt and the fudge factor amplitude Co = 10. 
The latter was obtained from the latest curve in Fig. |t] 
requiring that the value should be approximately 9 at 
x = 40mm. This is, of course, a rough and arbitrary 
choice regarding the strong noise in the plot, and its pur- 
pose is only to demonstrate the difference in the profile 
due to this procedure. 

It is obvious from Fig. llO] (b) that even for the rela- 
tively small ratio C b (0,t b ): C q (0, t q ) = 1/30, the pro- 
file is shifted to higher values for larger x, and that the 
expression given in Eq. ( |l9| ) seems to be a good approx- 
imation to the experimental result. It might be argued 
that this measure violates the requirement that the mean 
squared displacement should scale proportional to t 2 / 3 . 
We plotted the mean squared displacement correspond- 
ing to the modified profile, Eq. (19). It is obvious that 
the slope is approximately 2/3, due to the small rela- 
tive contribution of the Brownian solution (note that the 
integral determining the second moment of the Brown- 
ian contribution was cut off at x = 80mm since it has 
a non- negligible contribution for larger x). We expect 
that the relative amplitude of the Brownian contribution 
decreases with increasing system size, i.e., the finite size 
effects causing the necessity of the Brownian correction 
should become smaller. 



V. NUMERICAL EVALUATION BASED ON 
THE ORDINARY DIFFUSION EQUATION 

In addition to the analytical propagator treatment us- 
ing the fractional diffusion equation which describes the 
anomalous diffusion as an effective result of the complex 
boundary conditions, we applied numerical finite volume 
methods (FVM) to solve the ordinary diffusion equa- 
tion (^) for the explicit boundary conditions imposed by 
the geometrical structure of the finite model percolation 
cluster. The commercial software package FLUENT 5.5 
(TM) provides the numerical basis for this sort of analy- 
sis. 

The proton spin density distributions displayed in 
Fig. H were considered as a solution of the ordinary diffu- 
sion equation for the pore space of the percolation object 
and the initial condition of the experiment. Each lat- 
tice site is covered by a grid of 5 x 5 numerical unit 
cells. The diffusion process occurring with the ordi- 
nary diffusion constant of water at room temperature, 
D = 1.8 x 10~ 9 to 2 /s, was treated in a series of 20 s in- 
tervals. The proton concentration at the interface to the 



7 



reservoir was considered to be constant, C(0, t) — 1, for 
all times. 

Figure fl2| shows the normalized second moment (x 2 (t)} 
evaluated from the spatial proton density distribution as 
a function of time according to Eq. (12). As a conse- 
quence of the tortuous diffusion pathways in the percola- 
tion cluster, the time dependence turns out to be anoma- 
lous as expected. That is 



(x 2 (t)) oc *' 



0.8 



(20) 



The exponent value is coincides with the experimental 
result within the experimental accuracy, but is slightly 
higher than the theoretical expectation for random-site 
percolation clusters according to the Alexander/ Orbach 
conjecture, Eq. (||). 

The explanation is that there is a finite contribution 
of undisturbed Brownian diffusion trajectories through 
the pore space. The channel width relative to the exten- 
sion of the model system considered is not negligible as 
anticipated in the theoretical percolation cluster model. 
Therefore, a finite fraction of diffusion trajectories unaf- 
fected by the pore space restrictions contribute, and the 
exponent will be correspondingly larger. Of course, this 
"normal" contribution is expected to vanish in the limit 
of infinite system sizes. 



around 20 lattice constants at the occupation probabil- 
ities considered in this study (see Fig. ||). That is, the 
normal diffusion limit in Eq. (0) applies only extremely 
far above tj. This findingis in accordance with recent 
Monte Carlo simulations fig] , where the same conclusion 
was drawn. 

Apart from the mean squared displacement considered 
here, it is of interest to compare the whole concentra- 
tion profile with theoretical predictions. The concentra- 
tion profile in principle contains all information of the 
(anomalous) propagator effectively determining the dif- 
fusion properties. This is shown with the analytical prop- 
agator treatment based on the fractional diffusion equa- 
tion in Section IV. The experimental concentration pro- 
files for water diffusion in the model percolation clusters 
can be described very well on this basis. In this way, 
the fractional diffusion equation formalism was verified 
experimentally for the first time. 

Additionally, we have evaluated the ordinary diffusion 
equation for the boundary conditions given by the per- 
colation model cluster. The experimental data and the 
propagator results obtained from the fractional diffusion 
equation coincide completely as fas as can be judged in 
the frame of the accuracy of the diverse methods. 



VI. CONCLUSIONS AND DISCUSSION 
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FIG. 1: Photograph of a section of a quasi two-dimensional 
random site percolation model object (topview) (a) and of 
an entire model object (cross section) (b). The model object 
consists of several identical quasi two-dimensional percolation 
clusters stacked on each other in order to improve the signal 
intensity. The mechanical resolution of the fabrication pro- 
cess was 400 pm. The adjusted milling depth was constant 
between 1 to 2 mm in the various objects produced. 



FIG. 2: Proton spin density maps recorded in a two- 
compartment sample initially filled with bulk H2O and D2O 
gels (see the inset picture in Fig. []). The echo time was 
Te = 28 ms, the repetition time was Tr = 2 s. The field 
of view in x direction was 15 cm. The digital resolution was 
290 /Urn. The times indicate the span after contacting the two 
gels. The white lines represent the isotopic interdiffusion pro- 
files in the form of the projection of the proton spin density 
on the x direction. 



FIG. 3: Comparison of theoretical and experimental interdif- 
fusion profiles between bulk H2O and D2O gels. The data 
correspond to the two-dimensional spin density maps shown 
in Fig. m The solid lines represent fits of Eq. (|^) to the ex- 
perimental data. The fitted diffusion coefficients are given in 
the inset table. 



FIG. 4: Square of the position of the proton concentration 
profile at half height in bulk versus diffusion time. The exper- 
imental set-up consisted of two compartments initially filled 
with bulk H2O and D2O gels as shown in the inset picture. 
The largely asymmetric sizes of the H2O and D2O reservoirs 
ensure a practically constant proton concentration at the en- 
trance of the D2O compartment. The cross section of the D2O 
compartment was 1.5 cmx2 cm. The concentration profiles 
were obtained as projections of two-dimensional spin density 
maps on the x axis. The spin echo time was Te = 32 ms, the 
repetition time was Tr = 0.7 s. The field of view in x direc- 
tion was 7 cm. A digital resolution of 270 (jm was adjusted. 
The total acquisition time for a spin density map was 1 h. 



FIG. 6: Square of the position of the proton concentration 
profile at half height in the percolation cluster shown and de- 
scribed in Fig. ^ versus diffusion time. The inset shows the 
experimental set-up, where the percolation cluster (white) ini- 
tially is filled with D2O. The largely asymmetric sizes of the 
H2O and D2O reservoirs ensure a practically constant proton 
concentration at the entrance of the D2O compartment, i.e. 
the percolation cluster. The concentration profiles were ob- 
tained as projections of two-dimensional spin density maps on 
the x axis (see Fig. ||). The time resolution is 26'. The solid 
lines represent a fit of Eq. ^ for t <C t%. The fitted exponent 
parameter is d m — 2.89. 

FIG. 7: Typical isotopic inter-diffusion profiles directly mea- 
sured with the aid of one-dimensional imaging along the main 
diffusion direction. The experimental set-up is schematically 
shown in the inset of Fig. pi The spin-lattice relaxation time 
T\(x) was measured as a function of x by varying the rep- 
etition time Tr from 0.2 to 12.2 s in 15 stps. The signal 
intensity profiles were then corrected by multiplication with 
a factor (1 — exp{Tfl/Ti(a;)}) _1 , so that any spatially depen- 
dent saturation effects were eliminated. The spin-echo time 
was 20 ms, the field of view along the x direction was 15 cm. A 
diffusion time resolution of 1 h, and a digital space resolution 
of 290 /im were adjusted. 

FIG. 8: Mean squared proton displacement as a function of 
diffusion time in a two-dimensional random site percolation 
object the template of which is shown in the inset. The 
percolation cluster is initially (t = 0) filled with D2O. The 
characteristic data of the percolation cluster are: matrix size 
200 x 200; occupation probability relative to the threshold 
value p — p c = 0.030; fractal dimension df — 1.87; poros- 
ity of the percolating cluster p — 0.4845. The mean squared 
displacement data were evaluated from proton concentration 
profiles recorded in the form of one-dimensional NMR im- 
ages along the x axis. The proton concentration profiles were 
used to calculate the mean squared proton displacement in 
the percolation cluster moiety. The solid line represents a 
fit of Eq. (^) for t <C tj. The fitted exponent parameter is 
d w = 2.86. 



FIG. 5: Evolution of two-dimensional proton spin density 
maps (256 x 256 pixels) in a random site percolation object. 
The experimental setup is schematically shown in the inset of 
Fig. [| The times indicate the spans after attaching the H2O 
gel compartment to the multi-stack percolation model object 
(matrix size 100 x 100; p — p c = 0.029; df = 1.87, porosity 
of the percolating cluster p = 0.5352 ). The time resolution 
given by the image acquisition time was 26'. The echo time 
was Te = 23 ms, the repetition time was Tr = 0.7 s. The dig- 
ital resolution of the maps is 230 pm. The white lines overlaid 
to the spin density maps represent the mean proton concen- 
tration profiles (projections of the diffusion front on the main 
diffusion direction). 



FIG. 9: Comparison between the anomalous concentration 
profile 6*2/3(2;, t), left, and its Brownian analogue Ci(x,t), 
for a diffusion time interval t = 180 h. We used D2/3 = 
1.3 x 10~ 8 m 2 /s 2/3 and D = 2 x 10~ 9 m 2 /s. The amplitude 
is Co = 300 in both cases. The Brownian profile shows the 
much more efficient spread of the tracer substance into the 
medium. 

FIG. 10: (a) Anomalous concentration profile for a — 2/3, 
and with D2/3 = 1.3 x 10 _8 m 2 /s 2//3 . The times increase away 
from the origin, and are taken to be 20 h, 60 h, 100 h, 140 h, 
and 180 h, corresponding to Fig. ^. The inset shows a zoom 
into the plot range [0,25]. The conversion between t and 
the rescaled time r is given in Tab. |. (b) Corrected anoma- 
lous concentration profile, Eq. (|is|), with the subdiffusive pa- 
rameters Co = 300 and D 2 /3 = 1.3 x 10~ 8 m 2 /s 2/3 , and the 
Brownian diffusion constant D — 2 x 10~ 9 m 2 /s. The fudge 
amplitude for the Brownian contribution is 10. 
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FIG. 11: Mean squared displacement corresponding to the 
modified profile, Eq. (|l9|), log 10 -log 10 scale (full line). The 
dashed line represents the mean squared displacement cor- 
responding to the fractional result, Fig. ^. Both lines are 
approximately of slope 2/3. 



FIG. 12: Time dependence of the normalized mean squared 
displacement evaluated from numerical solutions of the ordi- 
nary diffusion equation for the boundary conditions given by 
the percolation model clusters (see Fig. ^|). The numerical 
procedure is based on the finite volume method (FVM). The 
numerical transient time resolution is 20 s. The time evolu- 
tion is scanned in steps of At = 3600 s. The data can be 
described by a power law (x 2 (t)} oc t ' 8 . 



t [h] 


/ \ 3 / 2 r -. ^ (\ 'X~\ 

(Ays) t [10- 6 m 3 ] 


Dt [10~ 4 m 2 ] 


20 


0.11 


1.44 


60 


0.33 


4.32 


100 


0.56 


7.20 


140 


0.78 


10.10 


180 


1.00 


12.96 



TABLE I: Conversion between experimental time scale and 
rescaled "time" t — Dt for Brownian and anomalous case. 
We use D — 2 x 10" 9 m 2 /s and D 2/3 = 1.3 x 10~ 8 m 2 /s 2/3 . 
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